SUBMITTED TO IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS 



1 



Block-Relaxation Methods for 3D 
Constant-Coefficient Stencils 
on GPUs and Multicore CPUs 

Manuel Rodriguez Rodriguez, Bobby Philip, Zhen Wang, Mark Berrill 



(N 

o 

< 



Abstract — Block iterative methods are extremely important as smoothers for multigrid methods, as preconditioners for Krylov 
methods, and as solvers for diagonally dominant linear systems. Developing robust and efficient algorithms suitable for current 
and evolving GPU and multicore CPU systems is a significant challenge. We address this issue in the case of constant-coefficient 
stencils arising in the solution of elliptic partial differential equations on structured 3D uniform and adaptively refined grids. Robust, 
highly parallel implementations of block Jacobi and chaotic block Gauss-Seidel algorithms with exact inversion of the blocks are 
developed using different parallelization techniques. Experimental results for NVIDIA Fermi GPUs and AMD multicore systems 
are presented. 

Index Terms — Parallel algorithms, Graphics processors, Smoothing, Multigrid and multilevel methods, Multicore 
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1 Introduction 

ITERATIVE methods such as Jacobi, Gauss-Seidel, 
Successive Over-Relaxation (SOR), and their vari- 
ants are extremely important as smoothers for multi- 
grid methods (2j, preconditioners for Krylov 
methods |3|, and as solvers for diagonally dominant 
linear systems. Their widespread use in iterative solu- 
tion methods for linear systems has led to significant 
effort being devoted to optimizing their performance 
on modern GPU and multicore systems. Much of the 
attention has been focused on point-wise versions of 
these methods to exploit parallelism GJ. 

However, block versions of these methods, in par- 
ticular line and plane smoothers, are extremely im- 
portant as key components of robust geometric multi- 
grid methods |5|, |6] and multilevel iterative methods 
(such as the Fast Adaptive Composite-Grid (FAC) 
and Multi-Level Adaptive Technique (MLAT) meth- 
ods 1 7 1, 1 8 1) on adaptively refined grids. Early work 
by Shortley and Weller |9] and Parter (To) has shown 
that even for isotropic diffusion problems block Ja- 
cobi and Gauss-Seidel can have advantages. More 
recently, Philip and Chartier demonstrated how to 
automatically construct block iterative methods for 
fairly general linear systems based on algebraic mea- 
sures of coupling (TT) . Furthermore, block iterative 
methods have the potential to increase local computa- 
tion and decrease communication on next generation 

• Manuel Rodriguez Rodriguez (rodriguezrom@ornl.gov) and Zhen 
Wang (1uangz@0rnl.gov) are with the Oak Ridge Leadership Comput- 
ing Facility, Oak Ridge National Laboratory, Oak Ridge, TN, 37831 

• Bobby Philip (philipb@ornl.gov) and Mark Berrill (berrillma@ornl.gov) 
are with the Computational Engineering and Energy Sciences Group, 
Oak Ridge National Laboratory, Oak Ridge, TN, 37831 



parallel systems where communication increasingly 
dominates costs. 

In this context, recent work focused on GPU im- 
plementations includes the 2D block based smoother 
of Feng et al. |12| and the ID block-asynchronous 
smoother for multigrid methods by Anzt et al. (13) . 
However, both works use Jacobi iterations within 
each block, rendering them closer to two-level Jacobi 
methods or Jacobi-Gauss-Seidel methods. These can 
be considered as variants on the work presented in 
Venkatasubramanian et al. [ 14 1 . Recent work in the 
context of multicore CPUs includes that by Adams et 
al. (151/ (16 1 on block Gauss-Seidel algorithms. 

In this paper we focus on developing efficient block- 
iterative Jacobi and chaotic Gauss-Seidel methods on 
current and evolving GPU and multicore architec- 
tures. We will demonstrate simple, efficient and gen- 
eral block smoothing algorithms with exact inversion 
of the blocks for 3D constant-coefficient elliptic prob- 
lems. Such problems are of importance in a wide 
variety of scientific computations. 

1.1 Model Problem 

We will focus on the model problem for the 3D 
Poisson equation 



V 2 w(x) = /(x), x= (x,y,z) e n, 

u (x) = 0, x e on, 



(i) 



for simplicity, though the presented methods will be 
applicable to any constant-coefficient elliptic system 
such as the recent work by Guy et al. |17] on block 
smoothers for Stokes problems. Here V 2 is the Lapla- 
cian operator, / is a source and u is the solution to 
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Fig. 1: In Figure 

from the shaded cells in Figure lb A domain composed of 4 patches is shown in Figure Id 



a patch is shown with an extraction of cells (shaded). The 7-point stencil in Figure 



lc 



arises 



Equation |lj on a cubic domain £1 = [0, l] 3 C M 3 with 
Dirichlet boundary conditions on the boundary d£l. 

It is assumed that the reader is familiar with dis- 
cretization methods and only a brief description is 
provided to establish notation. The methods pre- 
sented are not tied to any particular discretization, 
though for concreteness we use a cell-centered finite 
volume (FVM) discretization with variable unknowns 
located at cell centers. FVM discretization methods 
for single logically rectangular domains begin by 
partitioning the continuous domain fl into a set of 
discrete cells that together form a regular patch with 
n x , n y , and n z cells in the x-, y-, and z-directions 
respectively as in Figure [la] Each cell volume is then 
uniquely indexed by a triple (i,j,k) specifying loca- 
tion in space: [i ■ h x , (i + 1) • h x ] x [j ■ h y , (j + 1) • h y ] x 
[k ■ h z , (k + 1) ■ h z ] with < % < n x - 1, < j < n y — l, 
< k < n z — 1. Here h x ,h y , and h z represent the 
cell widths in each direction. To simplify notation 
going forward we assume h x = h y = h z = h. A cell- 
centered finite-volume approximation of equation |[TJ 
then leads to the following system of equations for 
each cell: 

6 u i,j,k — u i-l,j.k — u i+l,j,k ~ M»,j-l,k — Ui,j + l,k 

—lk,j,k-i — u ij,k+i — h 2 fij t k, (2) 

where Uij^ represents an approximation to u in cell 
(i,j,k). Since the coefficients in Equation |2j have 
no spatial dependence a stencil representation as de- 
picted in Figure lc can be used to represent both the 
coefficient and connectivity information for each cell 
in a patch. Cells on patch boundaries have the same 
stencil as interior cells, since we assume that each 
patch has a layer of 'ghost' cells around it (Figure |2b| 
to which boundary conditions are extrapolated. The 
methodology outlined extends with minor modifica- 
tions to discretizing a collection of non-overlapping 
logically rectangular patch domains as in Figure |ld| 
Interior ghost cells at patch interfaces are then used 
to ensure consistency across patches. 

Many simulations of physical phenomena often 
need to resolve extremely fine-scale features localized 
in space and/or time |18], |19|. Discretizing the entire 



(a) 



(b) 



Fig. 2: (a) Left: 3D AMR patch hierarchy with outlines 
of 8 patches on level and 8 patches on level 1 . Right: 
patch with shaded cells corresponding to a block, (b) 
2D patch slice with ghost and neighbor cells for a 
block (shaded). 



is then often impractical and not necessary. Instead 
Adaptive Mesh Refinement (AMR) can be used to 
increase the local resolution only where required. 
Given a collection of patches that cover the domain 
f^i = fl at some coarse resolution hi = h, structured 
AMR techniques identify a local subdomain 2 where 
finer resolution is required (f^2 may consist of dis- 
joint subdomains). A collection of patches with the 
same finer resolution h% are then introduced to cover 
f?2- Together, they form a refinement level covering 
the subdomain £1 2 - A typical choice of h 2 is hi/ 2. 
Repeating this process leads to a set of increasingly 
finer nested refinement levels, each consisting of a 
collection of logically rectangular patches. Patches in 
this hierarchy are dynamically created and destroyed 
as the simulation progresses depending on local res- 
olution requirements. Large-scale parallel structured 
AMR calculations can employ thousands of patches 
on each refinement level, with each processor owning 
multiple patches from potentially different refinement 
levels. Methods such as FAC and ML AT (8) smooth 
on these patches during the solution process. For 
the purposes of block smoothing, patches in the 
refinement hierarchy are typically too large to be 
processed efficiently. Instead patches are decomposed 
into smaller geometric blocks. For example, Figure 



2a 



physical domain at the required fine-scale resolution 



shows a 2 x 2 x 2 block (shaded) consisting of 
8 cells in a 4 x 4 x 4-cell patch. Thus, the patch of 
size 4 3 cells in Figure 2a is decomposed into 8 blocks 
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of size 2 3 . Stencil operations for each block in the 
patch have dependencies on neighbor and ghost cell 
data as shown in Figure [2b] As mentioned, the block 
relaxation methods we describe below are meant to 
be used as components of multigrid or multilevel 
methods in the two scenarios outlined of multiple 
structured patches over either a domain or a hierarchy 
of sub-domains in the case of AMR. 

The remainder of the paper is organized as fol- 
lows. Section |2] describes block-relaxation methods 
and algorithm choices based on GPU and multi- 
core CPU architectures. Section [3] discusses different 
parallelization strategies and considerations in the 
choice of block sizes. We state experimental results 
using state-of-the-art hardware including the NVIDIA 
TESLA X2090 GPU and a CRAY XE6 compute node 
in Section |4] The last section presents conclusions and 
an overview of possible future research directions. 

2 Block Smoothing Algorithms 

Discretizing equation ([TJ on a patch or a collection 
of patches leads to a linear system of equations, one 
equation per grid cell, of the form of equation Q. The 
resulting linear system of equations can be written in 
matrix form as 



Au 



(3) 



with A G W nxm , and u,f G W m with a suitable 
mapping from the matrix ordering to the (i,j,k) 
ordering on patches. Let A be partitioned into a set of 
submatrices as shown below. 



A = 



An 
A s i 



A 12 
Aii 



••• A u 

■ ■ ■ A 2s 

A -, A 

-^s.s — l -^ss 



(4) 



where Ay, 1 < i, j < s are now s 2 matrix subblocks 
of size qi x qj with J^i Qi — m - The partitioning of A 
into subblocks could be based on geometric blocking 
of a structured grid as described earlier or on other 
considerations such as anisotropic features in the PDE, 
or algebraic strength of coupling measures between 
variables. It is assumed that the diagonal blocks, 
An, 1 < i < s, are invertible. For the purposes of 
this paper it is sufficient to consider a lexicographical 
ordering of blocks within each patch. Furthermore, 
from this point on we will assume that all blocks are 
of the same size. A block stationary iterative method 
can now be defined by a splitting, A = M — N , where 
M is invertible, which leads to the stationary iteration 

u fe+1 = M~ l Nw k + M _1 f, (5) 

where k is the iteration number. The iteration given 
above converges if and only if p(M~ 1 N) < 1, where 
p denotes the spectral radius operator. An equivalent 
formulation of |5) is: 

„fe+i _ , 



for residual r k = f — Au fc , which will be the form used 
in this paper. An example of a block iteration is the 
block Gauss-Seidel method defined by the splitting 



M = 



An 

Asi 





A 22 



-4., 






As, 



(7) 



M" 



-i r fc 



(6) 



and N — M — A. When s = m, the subblocks are 
of size one and the iteration reduces to the standard 
lexicographic Gauss-Seidel iteration. A block Jacobi 
algorithm is obtained if Ay = for i > j also in 
equation |7|. The numerical results presented in this 
paper were all performed with block (chaotic) Gauss- 
Seidel and damped block Jacobi iterations, though 
once the block partitioning is defined we are free to 
choose any suitable block iterative process. 

In general, it is difficult to provide theoretical 
results that guarantee that a given block-iterative 
method formed will converge or that a particular 
choice of blocks will lead to a faster rate of conver- 
gence as opposed to an alternative choice of blocks. 
The interested reader is referred to Varga |20| for 
theoretical results on general block iterative methods 
and to Parter |T0) for some early work on block 
iterative methods for elliptic equations. 

The block-smoothing iteration derived from Equa- 
tion Q employs a direct solve to compute the inverse 
Ab := A^ 1 for the diagonal block A^ :— An with 
b,i = 1, . . . , s. The cost of computing this inverse for 
a block, and restrictions that apply in the case that 
the patch size is not a multiple of the block size, are 
described in Section 14.31 There it will be shown that 
it is sufficient to assume that only one inversion is 
necessary for a patch and this inversion has already 
been performed before the smoothing iterations are 
executed. 

Updating the neighbor cells of a block as in Figure 
2b can be performed in different ways. Our imple- 
mentation of the block smoothing algorithm employs 
either a classic Jacobi update or a chaotic block Gauss- 
Seidel update scheme between the blocks. This chaotic 
relaxation was first proposed by Chazan and Miranker 
|21| . Some current research on chaotic relaxation us- 
ing shared and distributed memory systems can be 
found under asynchronous iteration in p2j , |23) , pSJ . 
Baudet p3] defines the term chaotic relaxation scheme 
to describe a purely or totally asynchronous method, 
which accurately describes our chaotic block Gauss- 
Seidel implementation. Bertsekas et al. state in (25( a 
general convergence theorem for totally asynchronous 
algorithms in the case of fixed-point problems. A 
modified convergence theorem based on Chazan and 
Miranker with constraints on global memory and 
communication is presented by Strikwerda |26|. Blath- 
ras et al. state a timing model and stopping criteria for 
block asynchronous iterative methods in ]27|. A basic 
outline of the chaotic block Gauss-Seidel relaxation 
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for multiple patches either on an AMR refinement 
level or a multigrid level is shown in Algorithm [TJ 
assuming that the inner most loop will be executed 
asynchronously. In the block Jacobi update scheme 

Algorithm 1 Chaotic block Gauss-Seidel relaxation 
algorithm. 

for each step k do 
for each patch p do 
for each block b do 

u k b +1 = u k b + A b (f b - A b u k b ) 
end for 
end for 
end for 



a second vector v stores the new values of a patch 
during the smoothing step. At the end of the iteration 
the solution vector u for a patch swaps with v. This 
can be done in a very efficient way by using different 
pointers to the corresponding memory areas so that 
no additional memory copy is necessary. The chaotic 
block Gauss-Seidel update scheme does not need an 
extra vector. Advancing the solution can be done 
immediately. 

No synchronization between the blocks is necessary 
while processing the blocks within a smoothing step 
when using block Jacobi iteration or chaotic block 
Gauss-Seidel relaxation. Therefore the blocks can be 
processed independently without any communica- 
tion. The ghost cells between patches are updated in a 
Jacobi fashion at the end of a complete smoothing step 
over all patches. To prevent confusion about blocks 
on the GPU and blocks in the physical domain, a 
block in the domain is called a spatial block and a 
block on the GPU is called a thread block. The structure 
of AMR leads to the possibility to parallelize only 
over patches {patch parallel) or only over blocks {block 
parallel) leading to single-level parallelism. Parallelizing 
over patches and then over blocks leads to two-level 
or nested parallelism. 

A GPU is well suited to exploit the fine-grained 
parallelism typically exposed by single-level methods. 
It is also possible to exploit the parallel architecture of 
a GPU using the two-level parallelism by processing 
several patches and blocks concurrently. For two- 
level parallelism the number of patches that can be 
processed in parallel depends on the number of blocks 
that can be scheduled and executed independently on 
the GPU. This then depends on the block size, the 
hardware architecture of the GPU, and the number of 
blocks that are executed per patch. 

On multicore CPUs scalable performance is typi- 
cally achieved using coarse-grained parallelism. The 
work per thread must be large enough to amortize 
the overhead of scheduling the work. The block- 
parallel algorithm in which the blocks are processed 
in parallel has a potential bottleneck if the work per 



thread is low for small block sizes. Using the patch 
parallel method increases the workload per thread 
so that this bottleneck is avoided. Furthermore, two- 
level parallelism can improve the performance on a 
CPU beyond ordinary single-level parallelization |28j , 
which is also demonstrated in Section [4] 

2.1 Implementation for CUDA 

The implementation of block-relaxation algorithms on 
GPUs is done using the NVIDIA Compute Unified 
Device Architecture (CUDA). A short introduction 
into the technology and programming of CUDA can 
be found in |29|, |30|. Threads in CUDA are grouped 
together into thread blocks and are executed in a 
Single Instruction Multiple Threads (SIMT) fashion of 
32 threads called a warp. Thread blocks are organized 
into grids and each CUDA kernel launches one grid 
of several thread blocks. 

The decomposition of a patch in spatial blocks 
shares similarities with the CUDA architecture. But 
mapping a spatial block to a thread block or a patch 
to a grid on the GPU is nontrivial in the case of block 
smoothing. The mapping has an impact on the level of 
parallelism that can be achieved using different spatial 
block and patch sizes. We suggest the following two 
strategies for processing patches on the GPU: 

Strategy 1: One thread block processes all spatial 
blocks within a patch as in Figure 3a The thread 
block sweeps through the patch. In this case a GPU 
grid maps different thread blocks to different patches. 
Blocks in different patches are processed in parallel 
resulting in two-level parallelism. 

Strategy 2: One patch is processed by several thread 
blocks on the GPU as shown in Figure [3b] Here 
single-level parallelism exists if the patch size is 
large enough that the maximum number of scheduled 
thread blocks on the GPU is reached. In this case 
only one patch can be processed at a time on the 
GPU. If the patch size is small enough that multiple 
patches can be processed in parallel, then two-level 
parallelism is possible. 

In this paper we explore Strategy 2 for the presented 
implementations of the block smoothing relaxation. 
Strategy 1 will be explored in future research projects. 
In addition to the patch mapping, the spatial blocks 
must be mapped to the thread blocks on the GPU. 
A common method for spatial blocking is overlap- 
ping axis-aligned 3D blocking |31| , in which a spatial 
block of dimension b x x b y x b z with neighbor cells 
1Z is loaded into on-chip memory, leading to some 
redundant memory transfers for the neighbor cells. 3D 
blocking can be optimized by a 2.5D sliced blocking 
technique such as |32) or 1311. In 2.5D slicing, a 
block is decomposed into slices of cells. These slices 
are processed sequentially in an alternative blocking 
that preserves some cells in the on-chip cache or 
even registers, reducing global memory access. 2.5D 
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(a) 



(b) 



(c) 



(d) 



(e) 



Fig. 3: Examples demonstrating how to map patches to a GPU grid and how to map a spatial block to a 
thread block. Figure 3a shows several patches that are mapped to one GPU grid and |3b] shows one patched 



mapped to one GPU grid. In Figure 3e the thread block (dashed) is larger than the spatial block (shaded), 
in [3c] a thread block is smaller than the spatial block and in 3d a thread block has the same size as a spatial 
block. 



blocking is outside the scope of this paper but will be 
explored in future work. We focus on 3D blocking of 
the domain. Methods of creating a mapping from the 
spatial block to the thread block are: 

A. Thread blocks are larger than spatial blocks. Sev- 
eral spatial blocks can be processed by one thread 
block (Figure [3e|. 

B. Thread blocks are smaller than spatial blocks. The 
thread blocks have to sweep through each spatial 
block (Figure [3c|. 

C. Thread blocks are the same size as spatial blocks. 
One thread block processes one block in the 
spatial domain (Figure 3d I. 

The 2.5D slicing is the (B) block mapping, where the 
thread block sweeps through slices of a spatial block. 
Mapping (C) is used for the implementation of the 
block smoothing algorithm described in this paper. 
This mapping differs in two major ways from the 
other mappings. First, it is not necessary to sweep 
through the spatial domain, avoiding the introduc- 
tion of additional for-loops as would be necessary for 
mapping (B). Second, there is no need for a more 
fine-grained partitioning like in mapping (A). A fine- 
grained partition means that in order to identify a 
corresponding spatial block for a thread in mapping 
(A) additional indexing must be performed. 

In mapping (C) each thread in a thread block b is 
assigned a 3D index k) corresponding directly to 
the cell with the cell-centered spatial index. A thread 
resides in two index spaces: the global 3D space that 
is the index in the patch and the local 3D space that 
is the index within a thread block. The global 3D 
index calculation in CUDA without ghost cells using 
standard lexicographical ordering can be performed 
in the following manner: 

i = blockldx.x x blockDim.x + threadldx.x 
j = blockldx.y x blockDim.y + threadldx.y (8) 
k — blockldx.z x blockDim.z + threadldx.z. 

For the calculations with respect to the ghost cells, the 
thread index must be shifted by the ghost cell width. 



The local 3D index is given by the three threadldx 
values in CUDA. 

Cells of a patch including the ghost values are 
stored in a ID array. The right-hand side of Equation 
|[3} needs no ghost values. Accessing both vectors on 
the GPU kernel requires different indexing. Thus a 
thread resides in additional ID index spaces: the global 
ID space with and without ghost values. If the on-chip 
shared memory is used within a block, an additional 
local ID space is also required. The index calculation 
for the ID space is a simple axis projection either on 
the local blockDim or the global blockDim x gridDim 
axis. In the resulting block Jacobi iteration on the 
GPU, as shown in Algorithm [2| each thread in a 
thread block b computes the new value of a cell in 
the physical domain. In the first step of Algorithm [2] 

Algorithm 2 Block Jacobi algorithm on GPUs. 

{Step 1 - Thread block computes residual} 

n =fb- A b u b 

{Barrier for shared memory update} 
syncthreads() 

{Step 2 - Thread block computes dgemv} 

r b = A b r b 

{Step 3 - Thread block advances solution} 

v b — u b + ur b 

{Step 4 - Swap the vectors} 



the computation of the block residual r b is performed, 
including the 7-point stencil operation from Equation 
|2j that each thread performs. It is possible to store 
the constant stencil in constant memory on the GPU. 
The constant memory on the GPU is a fast, read-only 
cached memory area in the global memory space of 
the GPU. The computed residual r b for a particular 
spatial block b is stored in shared memory on the GPU 
because all threads in a block must access the residual 
in order to compute the matrix-vector operation in 
step 2. 

The second step of Algorithm [2] consists of multi- 
plying the block-diagonal inverse A b £ M. nxn (where 
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n is the product of the block dimensions) with the 
block residual r&. Each thread in a thread block per- 
forms a dot product for its particular cell. The matrix 
is stored in column-major ordering to optimize the 
memory access pattern. Step 3 advances the solution 
of a spatial block. In the case of Jacobi iteration the 
relaxation parameter lu = 0.8 is used and the new 
block solution is stored in a second vector v. The 
chaotic block Gauss-Seidel updates the solution to 
u b + ujr b with uj = 1 without any intermediate vector 
v. Swapping the vector v with u in step 4 is only 
necessary for the block Jacobi iteration. An implemen- 
tation of the chaotic block Gauss-Seidel CUD A kernel 
is reproduced in the Appendix. 

A feature of the NVIDIA GPUs is that the GPU sup- 
ports concurrent data transfers, transfer-execution or 
kernel-execution overlapping |33|. Depending on the 
instruction stream in the algorithm, either concurrent 
data transfers, transfer-execution or kernel-execution 
overlapping is performed. 

For transfer-execution overlapping, the instructions 
must be executed in depth-first order where for all 
patches the instruction chain: memory transfer to the 
GPU, kernel execution and the memory transfer back 
to the host is executed. In this case the transfers from 
and to the GPU memory can overlap with kernel exe- 
cution. To execute kernels concurrently, the sequence 
of instructions must be in breadth-first order in which 
the transfer for all patches is started, then the kernel 
execution for all patches and at last the transfer back 
to the host for all patches. The overlapping of kernel 
executions is possible only when the number of blocks 
executed by one kernel does not exceed the number 
of blocks that can be scheduled concurrently on the 
GPU. 

2.2 Implementation for a Multicore Architecture 

In the serial implementation of block Jacobi, Algo- 
rithm [3] will be executed for each patch. The CPU 
implementation first computes the block residual. Af- 
ter this step the block residual is multiplied with 
the block inverse in step 2. This is done using the 
optimized Basic Linear Algebra Subprograms (BLAS) 
Level-2 (dgemv) function. Step 3 must be executed 
separately from step 2 because accessing the vector 
v requires a global 3D index computation. Finally 
advancing the solution using the block Jacobi update 
scheme has to be done after all blocks are processed 
in step 4. 

Algorithm [3] is parallelized with OpenMP (34j di- 
rectives by implementing the single-level and two- 
level parallelism. OpenMP is a compiler extension 
that enables portable, automated thread parallelism 
for multicore CPUs. Additional information about 
OpenMP programming can be found in |35|, |36|. 
For the single-level parallelism the implementation 
parallelizes either the processing of patches or the 



processing of blocks. In the two-level parallelism p 
threads process the patches and spawn additional q 
threads to process the blocks. 

An OpenMP parallel for loop is used to parallelize in 
the single-level parallelism over either the blocks or 
the patches. For the two-level parallelism an OpenMP 
parallel region is created to parallelize over the patches 
and an additional region is created to parallelize over 
the blocks using the OpenMP 3.0 collapse directive to 
collapse the loops over the 3 dimensions x,y,z. 



Algorithm 3 Block Jacobi algorithm on CPUs. 

for block index in Z do 
for block index in Y do 
for block index in X do 

{Step 1 - Compute block residual} 

n = h - A b u b 

{Step 2 - Compute dgemv} 

r b = A b r b 

{Step 3 - Create block solution} 
v b = Lur b 
end for 
end for 
end for 

{Step 4 - Advance the patch solution} 

U = U + V 



The chaotic block Gauss-Seidel version of Algo- 
rithm [3] uses only the single-level block parallelism 
or nested parallelism to process the blocks asyn- 
chronously. 

3 Estimating Spatial Blocking 

The block size in the physical domain influences 
the convergence and smoothing properties of the 
algorithm. For multilevel algorithms it is the latter 
criteria that is of more importance and is also less well 
understood theoretically other than for specific cases. 
However, it is important that the algorithm converge 
all error modes and our numerical experiments in this 
subsection demonstrate this. 

The results of numerical experiments on conver- 
gence behavior with different block sizes for block 
Jacobi and chaotic block Gauss-Seidel relaxations are 
presented in Figure |4] The following block sizes 
(b x ,b v ,b z ) are used: (2,2,2) = 8, (4,2,2) = 16, 
(4,4,2) = 32, (4,4,4) = 64, (8,4,4) = 128, (8,8,4) = 
256, (8,8,8) = 512. We use (4,2,2) instead of (2,4,2) 
or (2,2,4) because the x-dimension provides the 
fastest memory access in our implementation. We note 
that for problems with anisotropic coefficients this 
may not be the best choice. 

For the convergence study a patch size of 64 3 is 
used. The relative error estimate is computed after 3 
smoothing steps. In multigrid applications, perform- 
ing at most 3 smoothing steps is usually sufficient to 
smooth the oscillatory error modes. The selection of 
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(a) Relative error. (b) Convergence factor. 

Fig. 4: Convergence study of block Jacobi and chaotic block Gauss-Seidel relaxation. 



the 64 3 patch size is based on the volume to surface- 
area ratio and the overall memory required to store 
a patch. Considering a ghost cell width of 1, a 32 3 
domain with ghost cells needs overall 34 3 cells. The 
memory overhead to store the ghosts cells is 6536 
cells, around 20% of the overall memory needed for 
32 3 cells. By using a patch size of 64 3 the overhead 
is only 10% and for 96 3 only 6%. Larger patches 
need more memory space so that a patch size of 96 3 
requires about 14 Megabytes in double precision to 
store u and / for all cells in the Gauss-Seidel algorithm 
and about 21 Megabytes in the Jacobi algorithm. With 
bigger patches the number of patches that can be 
processed within a GPU or CPU drops which impacts 
the balance of the patch distribution across different 
processes. Choosing patches between 64 3 and 96 3 is 
a trade-off between the number of patches, memory 
usage and load balancing. 
In Figure 



4a 



the relative error ||rfc|| / ||ro|| is shown 



for different block sizes for k = 3. Figure 4b shows 
the asymptotic behavior of the convergence factor 
v/||rj_|_fc|| / ||rj|| where i, k defines a sliding window 
for a random initial guess and a large number of 
smoothing steps (« 10000). The asymptotic conver- 
gence rate for block Jacobi is reached after a relatively 
few number of steps while in the case of chaotic 
block Gauss-Seidel the rate will jitter from step to step 
because of the chaotic updates but nevertheless the 
estimation is fairly accurate. 

The block Jacobi and chaotic block Gauss-Seidel 
convergence study shows an overall monotonic im- 
provement in relative error and convergence factor 
by increasing the block size for a random initial 
guess. The convergence in smoothing is necessary 
but does not reflect in general which types of modes 
have been smoothed and is outside the scope of this 



paper. However, the convergence behavior serves as a 
baseline for the performance evaluation with respect 
to the block size for the presented implementation. 

4 Experimental Results 

This section presents the experimental results for var- 
ious GPU and multicore CPU implementations of our 
block smoothing algorithm. In Section 4.1 the speedup 
of a 24-core system will be presented. The speedup 
S(p) of p processors is defined as the fraction of 
sequential over parallel wall time on p processors. 
Estimating S(p) with a constant problem size and 
increasing processor count is called strong scaling. 
A speedup of S(p) = p is ideal, as the wall time 
improves proportional to the number of processors 
that are added. Strong scaling gives a good under- 
standing of how well the algorithm scales with the 
tested parameter set. We can derive the efficiency 
E(p) of a parallel system with p processors from the 
speedup. The efficiency is defined by the ratio S(p)/p 
between the speedup and the number of processors 
that reaches the particular speedup. 



Section 4.2 compares the wall time of different pa- 
rameter sets executing the block-relaxation algorithms 
on the GPUs and CPUs. For the benchmarks different 
hardware including a CRAY XK6 node (TitanDev) 
from the "Jaguar" supercomputer at the Oak Pudge 



National Laboratory are used. Finally in Section 4.3 



we provide detailed benchmarks for inverting the 
diagonal-block matrices of equation ||7| with different 
block sizes. 

All measurements presented in this section do not 
include the update of the ghost cells between patches 
because our focus is on the stand-alone kernel execu- 
tion time. However, the wall time includes the transfer 
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Fig. 5: Speedup of block Jacobi and chaotic block Gauss-Seidel relaxation using different levels of parallelism. 



from and to the GPU for all patches. We observed less 
than 10% of the best wall time on the GPU and around 
5% of the wall time on the CPU for updating the ghost 
cells. This time grows linear with the number of cells. 

4.1 Multicore Speedup 

A state-of-the-art Symmetric Multiprocessing (SMP) 
compute node equipped with double AMD Opteron 
6128 Magny Cours CPUs is used to benchmark the 
multicore implementations. Each CPU contains 2 dies 
and each die has 6 floating-point units. Using the 
two-level parallelism 4 threads (number of dies that 
the system has) are spawned to process the patches 
in parallel and overall 24 threads are processing the 
blocks in parallel. Figure 5] shows the strong-scaled 
speedup for the block Jacobi and chaotic block Gauss- 
Seidel algorithm processing 96 patches with patch size 
of 64 3 using two different block sizes. The selected 



block sizes are 4 3 and 8 3 because the experiments in 
Section [3] suggest a good convergence factor for 8 3 . 
The block size of 4 3 is used to provide a contrast to 
the 8 3 with respect to the cost of the matrix-vector 
multiplication, the inversion of the block diagonal 
matrix and performance on the GPU and CPU. 

The results from the multicore benchmarks show 
that using either the patch-parallel or two-level par- 
allelism gives the best scalability for all cases. A 
maximum speedup around 18 is reached for the block 
Jacobi and chaotic block Gauss-Seidel using a block 
size of 8 3 with two-level parallelism. 

This leads to an efficiency of about 75% on the 
24-core machine. Choosing a block size of 4 3 drops 
the maximal speedup for both algorithms to about 
12 for block Jacobi and 14 for chaotic block Gauss- 
Seidel. The patch-parallel version of the chaotic block 
Gauss-Seidel is in fact a standard block Gauss-Seidel 
relaxation without chaotic updates because the blocks 
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Name 


Processors 


Cores 


Speed 


GFLOPS 


Watts 


Interlagos (XE6) 


2 x AMD 6272 


4 x 8 = 32 


2.1 GHz 


-295 


230 


Maranello 


2 x AMD 6168 


4 x 6 = 24 


1.9 GHz 


-202 


230 


Fermi 


1 x NVIDIA C2050 


14 x 32 = 448 


1.2 GHz 


-515 


238 


Fermi+ (XK6) 


1 x NVIDIA X2090 


16 x 32 = 512 


1.3 GHz 


-655 


225 



TABLE 1: Test systems for the performance benchmarks. The column Cores shows the processor core 
configuration, the columns Speed and GFLOPS the chip clock speed and the theoretical peak performance 
in double precision of the cores and the column Watts the power consumption in Watts. 
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Fig. 6: Chaotic Gauss-Seidel relaxation on Fermi+ using different overlapping techniques. 



are processed in serial. 

A larger block size leads to better scalability with 
respect to efficiency for both algorithms. The two-level 
parallelism scales as well as the patch-parallel version 
and allows the chaotic block update in the Gauss- 
Seidel scheme. 

From the speedup results we see that the block 
Jacobi algorithm gains more from processing the 
patches in parallel than processing the blocks in paral- 
lel, because the block-parallel version does not scale 
well with smaller block sizes. The patch-parallel or 
nested-parallel versions of Jacobi scale better because 
the solution u k+1 is advanced for the patches in par- 
allel. In the case of the chaotic block Gauss-Seidel the 
block parallel version scales better for smaller block 
sizes than the two-level parallelism. The stairstep 
pattern in the speedup of the two-level parallelism 
indicates a workload imbalance caused by adding 
block parallelism to the patch parallel version. 

4.2 Algorithm Performance 

In this section we present measurements that have 
been taken on the systems listed in Table 1 and from 
now we reference the system by its name. One system 
named 'Maranello' consists of two AMD Opteron 
6128 CPUs and one NVIDIA C2050 (Fermi) (37). The 



other system named 'Interlagos' is a CRAY XE6 |38) 
node equipped with two AMD Opteron 6272 CPUs 
and a CRAY XK6 node 139] provides the bench- 
marks for the NVIDIA TESLA X2090 (Fermi+) GPU. 
Floating-Point Units (FPUs) on the NVIDIA Fermi 
GPUs are arranged as a set of 32 in each Streaming- 
Multiprocessor (SM). The major differences between 
the two GPUs are that the Fermi has 14 SMs and 
the Fermi+ is equipped with 16 SMs and the Fermi+ 
GPU has a higher GPU and memory clock speed. On 
the CPUs the arrangements of FPUs are different as 
for Maranello each die consists of 6 FPUs and for 
Interlagos each die has 4 FPUs but 8 integer units. 
The different FPU arrangements reflect one difference 
between the two AMD architectures. All benchmarks 
in this section exclude the time to compute the in- 
verses. Timings to compute the inverses are presented 
in Section 14.31 



The GPU has the option to overlap as described 
in Section 2.1 A benchmark is used to explore the 
performance of the different overlapping techniques 
to estimate which method should achieve the best 
performance for different patch sizes. Figure [6] shows 
the wall time for 3 smoothing steps varying numbers 
of patches and the patch size on the Fermi+ by 
exploiting kernel overlapping (breadth-first) or kernel 
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Fig. 7: Time to smooth with chaotic block Gauss-Seidel and block Jacobi relaxations. The data is grouped by 
parameter set (patch size / number of patches / block size) for different systems. 



transfer overlapping (depth-first). Figure [6] also shows 
the wall time for the Interlagos system using the two- 
level parallelism for comparison. 

One patch of size 16 3 does not fully occupy the 
GPU, so multiple patches can be smoothed concur- 
rently (two-level parallelism) using the breadth-first 
approach. For a 16 3 patch size with a block size of 
8 3 a total of 8 blocks are executed. A thread block 
is scheduled on a SM, and 8 thread blocks can be 
active at a time on a Fermi SM so that only one SM is 
occupied by a patch, leading to 16 concurrent patches 
on a Fermi+ and 14 on a Fermi GPU. Using breadth- 
first approach with smaller patches is faster than 
depth-first as stated in Figure 6a With kernel transfer 
overlapping (depth-first), the GPU does not process 
the next patch until the current patch is finished, 
leading worse runtimes for smaller patches. 

Increasing the patch size changes this behavior as 
shown in Figure 6b For a patch of size 64 3 and a block 
size of 8 3 , a total of 262,144 threads are spawned in 
512 thread blocks so that each SM on a Fermi+ is over- 
subscribed by a factor of 4 with this configuration. 
Therefore, no concurrent kernel execution is possible 



instead overlapping kernel execution with memory 
transfer is performed. Thus, the memory transfers for 
the patches can be hidden completely by using the 
depth-first approach with a patch size of 64 3 , giving 
results superior to the breadth-first approach. 

Decreasing the block size will not improve results, 
because the total number of active threads on the GPU 
is decreasing. For example by using a block size of 4 3 , 
only 512 threads are active on a SM. On a Fermi SM 
48x32 = 1536 (Active warps x Warp size) threads can 
be active at a time so that a block size of 4 3 reaches 
only 1/3 of the maximum number of threads. The 
limiting factor in this case is the block size. Choosing 
a block size of 8 3 saturates a SM just by executing 3 
blocks, reaching over 90% of measured occupancy. For 
the NVIDIA Kepler |40| architecture the maximum 
number of active threads is 2048, so 4 blocks can be 
active at a time on a Kepler Streaming-Multiprocessor. 

For the benchmarks the focus is on patch sizes 
of at least 64 3 so that the timings are taken with 
depth-first asynchronous memory transfer from and 
to the GPU in each smoothing step, while the CPU 
implementation does not need to execute memory 
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Fig. 8: Smoothing sets of patches with the same and with mixed sizes. 



transfers after each step. Timings for the CPU code 
are taken from the two-level nested version. In Figure 
[7] the wall time, smoothing selected parameter sets 
is shown for the block Jacobi and for chaotic block 
Gauss-Seidel relaxations. 

Three smoothing steps are executed on varying 
parameter sets. Each parameter set is defined by patch 
size, number of patches and block size. Sets with 
patch sizes of 64 3 and 96 3 and block sizes of 4 3 and 
8 3 are benchmarked. The Interlagos outperforms the 
Maranello because of its higher peak performance. 
Both algorithms scale with an efficiency near 75% on 
the Interlagos, only if 16 threads are used and the 
threads are placed on alternating cores to avoid multi- 
ple threads sharing one floating point unit. Beyond 16 
threads the efficiency drops to 50-60% but the overall 
wall time still improves. 

Both GPUs perform better than the multicore sys- 
tems. The Fermi+ is faster then the Fermi because of 
the higher peak performance and the additional cores 
on the Fermi+. The wall time using the parameter set 
(64 3 /96/8 3 ) for both algorithms on the Fermi+ is about 
1.8 x faster than the multicore time on the 32 cores per 
node CRAY XE6 (Interlagos) machine. 

The selected application parameter sets use the 
same patch sizes which are characteristic of one class 
of AMR methods. However, for another class of AMR 
applications the patch size can vary We benchmark 
the performance of the block smoothers using dif- 
ferent patch sizes to better reflect the practical per- 
formance for this latter class of AMR methods. For 
estimating the practical performance, a benchmark 
executing the smoothing algorithms on patch sets 
with the same patch sizes and sets with different sizes 
as listed in Table 2. All sets consist of 80 patches 
and the mixed patch set is divided into subsets of 



16 patches with 64 3 , 72 3 , 80 3 , 88 3 and 96 3 . The ratio 
of cells between 64 3 and 80 3 is 1.89 and the ratio of 
cells between 64 3 and the mixed set is 2.07. The ratio 
of cells between 80 3 to 96 3 is 1.72 and the ratio of cells 
between mixed an 96 3 is 1.63. This ratio is important 
because there exists a correlation between the number 
of cells and the wall time for smoothing the cells. 



Patch size 


#Patches 


#Cells 


Block size 


64 x 64 x 64 


80 


62,914,560 


8x8x8 


80 x 80 x 80 


80 


122,800,000 


8x8x8 


Mixed 


80 


130,252,800 


8x8x8 


96 x 96 x 96 


80 


212,336,640 


8x8x8 



TABLE 2: Patch sets for smoothing fixed number of 
patches using sets of the same and mixed sizes. 

Figure [8] shows the wall time for smoothing patches 
of the same and mixed sizes using block Jacobi or 
chaotic block Gauss-Seidel on the fastest GPU and 
the CPU system. The time for smoothing the patches 
on the GPU increases roughly proportionally to the 
number of cells independent of the mixed patch sizes 
for both algorithms. This behavior is expected be- 
cause no synchronization is needed at the end of a 
processed patch and the memory transfer overlaps 
with the kernel execution. Essentially the GPU just 
batch processes block after block without stalling. 
However, the behavior is not the same in the CPU 
implementation. In the CPU implementation the wall 
time increases proportionally to the patch size for 
constant-size patches. For the mixed patch size set 
the wall time is closer to the time that is needed to 
smooth the largest patches. The reason for the poor 
performance of the mixed set is that the OpenMP 
threads looping over the blocks must join after the 
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patches are processed. 

Load balancing in the case of mixed patch sizes is 
an issue that depends on the ordering of the patches. 
In reality, the ordering of patches is random and in 
the worst case one thread gets all the larger patches. 
Using a naive approach by partitioning the patches 
without respect to the size has a potential for imbal- 
ance. The computation of an optimal partitioning is a 
variation of the bin packing problem which is known 
to be NP-hard. We can control the ordering of patches 
in the benchmark and if the patches are partitioned 
perfectly (each thread has the same amount of work) 
then the wall time for the mixed set drops to 3.576 
seconds for block Jacobi and 3.409 seconds for chaotic 
block Gauss-Seidel, leading to a wall time to number 
of cells relation like the GPU. 

Smoothing using mixed patch sizes accelerates the 
GPU wall time more than 1.8 x over the CPU wall 
time in the case of block Jacobi relaxation and chaotic 
block Gauss-Seidel when the patches are imbalanced. 
A balanced patch partitioning improves the wall time 
on the CPU for the mixed set, but to be fair the 
additional cost for partitioning must be taken into 
consideration because it is not necessary on the GPU. 



O 
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a> 
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Block size 

Fig. 9: Kernel performance for chaotic block Gauss- 
Seidel using different block sizes with patch size 64 3 . 

The correlation between wall time and patch size 
or wall time and number of patches can be clearly 
derived from the stated benchmarks in Figure [8] and 
Figure [6j But as we showed in Section [3| another 
important factor contributing to the wall time is the 
block size. The number of floating-point operations 
for the stencil is constant for a thread in a thread 
block regardless to the block size, so the performance 
is dominated by the matrix-vector multiplication of 
the block inverse with the residual. Figure [9] shows 
the wall time for smoothing 3 iterations on 64 patches 
with a patch size of 64 3 using different block sizes for 
the chaotic block Gauss-Seidel method. 

The CPU performs better for the block size of 2 3 . 



As the block size increases, the GPU outperforms the 
CPU. For a block size of 8 3 more than 56 million 
cells per second can be smoothed on the Fermi+ and 
around 28 million cells per second on the Interlagos. 

4.3 Block inversion 

Our block-relaxation algorithms require block inver- 
sions of diagonal blocks from the matrix that corre- 
sponds to the AMR domain. In the case of constant- 
coefficient stencils the computation of the inverse 
must be performed only once if the patch size is a 
multiple of the block size. If the patch size is not a 
multiple of the block size then two possibilities exist. 
Either the blocks that do not fit at the end of the 
patch overlap with the previous ones or all inverses 
for these cases must be precomputed. 




32 64 128 

Block size 

Fig. 10: Total time to compute the inverse of either one 
matrix or all block matrices for different block sizes. 

Overlapping can potentially change the conver- 
gence behavior, even resulting in divergence of the 
smoothing algorithm for chaotic updates. Computing 
all inverses is a robust approach that needs more 
memory and more computation time upfront. A total 
of b n = b x b y b z inverses are potentially needed for 
a block of size b x x b y x b z . To estimate the inverse 
computation overhead a benchmark computes either 
one inverse for a block or all inverses for a particular 



block size as shown in Figure 10 



The inversion is performed using a sequential BLAS 
library on the Interlagos system. Computing the in- 
verses can be done in parallel but the cost of com- 
puting the inverses in serial is negligible in both 
cases. For example, computing the inverse of the block 
size 8 3 takes 0.06 seconds on Maranello, or about 
5% of the wall time of 1.19 seconds needed for one 
simulation iteration smoothing three steps with 64 
patches of size 64 3 using chaotic block Gauss-Seidel 
on the same system. In the case of multiple inverses 
this fraction increases but is also negligible because 
the inverses have to be computed only once for the 
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whole simulation in the case of constant-coefficient 
stencils and therefore the time for precomputing the 
inverses is amortized over a few simulation steps. For 
example, in our typical AMR simulations the inverses 
can be reused thousands of times. 

5 Conclusions and Future Work 

We implemented our block smoothing algorithm for 
structured adaptively refined meshes using block Ja- 
cobi and chaotic block Gauss-Seidel relaxations for 
modern GPUs and multicore CPUs. Our multicore 
implementation scales on a 24-core machine with an 
efficiency of 75%, achieving a speedup around 18 
by exploring multiple parallel strategies. The GPU 
version is about 1.8x faster than our best multicore 
system using patches of the same sizes. In practical 
benchmarks with mixed patch sizes we can observe 
the same wall time improvement using a GPU. The 
peak performance of the X2090 is about 2.2x the peak 
of the newest CRAY XE6 compute node. An overall 
wall time acceleration from 1.5x-2x on the X2090 
is a very good and realistic result but there is still 
room to improve both implementations. Furthermore, 
we expect an additional improvement of the GPU 
performance with the NVIDIA Kepler architecture. 

The way asynchronous memory transfer on the 
GPU will be issued has an impact on the wall time 
for small patches. Using depth first computation with 
transfer-overlapping gives the best result for large 
patches and can hide the memory transfer completely 
so that no disadvantage of the memory transfer occurs 
against the CPU version where no memory transfer is 
needed. Comparing the block Jacobi with the chaotic 
block Gauss-Seidel shows that both algorithms scale 
in the same fashion and the chaotic block Gauss- 
Seidel algorithm gives a better wall time. The chaotic 
block Gauss-Seidel gives better results in terms of 
convergence and relative error for the random initial 
guess across all tested block sizes. Increasing the block 
size further improves the convergence factor for both 
algorithms. 

In addition we presented strategies for the paral- 
lelization of the block smoothing algorithm in the con- 
text of adaptively refined meshes. Different abstract 
mapping strategies for the GPU were also presented 
and can be used independently from the smooth- 
ing in future structured AMR-based algorithms. The 
implemented OpenMP versions use single-level and 
two-level parallelism approaches to realize these map- 
pings. 

Our future investigation is to optimize our block- 
relaxation implementations for the next hybrid super- 
computer architecture by exploring more GPU par- 
allelism like tiling the blocks or other patch map- 
ping strategies. Moreover a hybrid implementation by 
smoothing patches on the GPU and CPU is necessary 
because the memory transfer and kernel execution on 



the GPU is asynchronous so that the CPU is idle dur- 
ing the processing on the GPU. This wasted compute 
power can be used in future implementations. 
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Appendix 

cuda source for chaotic block 
Gauss-Seidel using a 7-point stencil 



global void smooth3d7c (A, u, f) { 

extern shared double sres[]; 

// Global 3D index in the physical domain 

int i = blockldx.x * blockDim . x + threadldx.x; 

int j = blockldx.y * blockDim . y + threadldx.y; 

int k = blockldx.z * blockDim . z + threadldx.z; 

// Global ID index in the physical domain 

int gid = k * gridDim . x * blockDim . x 

* gridDim. y * blockDim. y 

+ j * gridDim. x * blockDim . x + i; 

// Local ID index in the thread block 

int lid = threadldx.z * blockDim . y * blockDim . x 

+ threadldx.y * blockDim . x 

+ threadldx.x; 

// Block size in device/physical domain 

int n = blockDim. x * blockDim. y * blockDim . z ; 

int nx = gridDim. x * blockDim . x + 2; 

int ny = gridDim . y * blockDim. y + 2; 

int xy = nx * ny; 

// Set residual with right hand side 
double r = f[gid]; 

// Shift index by ghost cell width of 1 
k += 1; j += 1; i += 1; 
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// Apply stencil 

r — = stencil [6] * u[(k+l)*xy + (j)*nx + (i)]; 
r — = stencil[5] * u[(k — l)*xy + (j)*nx + (i)]; 
r — = stencil[4] * u[(k)*xy + (j+l)*nx + (i)]; 
r — = stencil[3] * u[(k)*xy + (j — l)*nx + (i)]; 
r — = stencil [2] * u[(k)*xy + (j)*nx + (i+1)]; 
r — = stencil [1] * u[(k)*xy + (j)*nx + (i— 1)]; 
r — = stencil[0] * u[(k)*xy + (j)*nx + (i)]; 
gid = (k)*xy + (j)*nx + (i); 
sres [ lid ] = r ; r = 0; 
syncthreads ( ) ; 

// Multiply the block with the inverse 
#pragma unroll 
for(int g = 0; g < n; g++) 
r +=A[g*n+lid] * sres[g]; 

// Advance solution 

u[gid] = omega * r + u[gid]; 

} 
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